#include "topology.h"
#include "combination.hpp"
// Two-loop topologies t1 and t2

Graph t1_full()
{
  Graph t1;
  add_edge(0, 1, t1); 
  add_edge(1, 3, t1); 
  //  add_edge(1, 3, t1); 
  add_edge(3, 2, t1); 
  add_edge(2, 0, t1); 
  add_edge(2, 1, t1); 
  return t1;
}

 const char* topo_map[] = 
{
// one-loop topologies
  "l1",              // 
  //two-loop topologies
  "t1",            //
  "t2",            //
  //three-loop topologies
  "la",              //
  "be",
  "bu",
  "fa",
  "no",
  "y1",
  "y2",
  "y3",
  "y4",
  "y5",
  "o1",
  "o2",
  "o3",
  "o4"  
};

Graph_topo get_root_topos(int l)
{
  BOOST_ASSERT_MSG((1<=l)&&(l<=3),"only 1,2,3-loop topologies are defined");
  Graph_topo out_vec;
  bool inserted;
  Graph::edge_descriptor edis;
  switch(l)
    {
    case 1:
      {
//************* Topology l1***********
//               
//            ____/____
//           /    \    \   
//          /           \ 
//   [Q]___/0           1\___
//         \             /
//    [P4]  \           / 
//           \____\____/
//                / 
//
//***********************************
        Graph l1;
        tie(edis,inserted) = add_edge(0,1, l1); 
        l1[edis].gmomentum = get_symbol("p1");
        l1[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(1, 0, l1); 
        l1[edis].gmomentum = get_symbol("p2");
        l1[edis].gexpr = 1/get_symbol("(p2.p2)");
        out_vec.push_back(std::make_pair(l1,Tl1));
        break;
      }
    case 2:
      {
//************* Topology T1***********
//               
//            ____3____
//           /    |    \_  
//     [P1]|/     |    |\ [P2]
//   [Q]___/0 [P5]|     2\___
//         \     /|\    _/
//    [P4] _\|    |     /| [P3]
//           \____|____/
//                1
//
//***********************************
      Graph t1;
      tie(edis,inserted) = add_edge(3, 0, t1); 
      t1[edis].gmomentum = get_symbol("p1");
      t1[edis].gexpr = 1/get_symbol("(p1.p1)");
      tie(edis,inserted) =   add_edge(2, 3, t1); 
      t1[edis].gmomentum = get_symbol("p2");
      t1[edis].gexpr = 1/get_symbol("(p2.p2)");
      tie(edis,inserted) =   add_edge(1, 2, t1); 
      t1[edis].gmomentum = get_symbol("p3");
      t1[edis].gexpr = 1/get_symbol("(p3.p3)");
      tie(edis,inserted) =   add_edge(0, 1, t1); 
      t1[edis].gmomentum = get_symbol("p4");
      t1[edis].gexpr = 1/get_symbol("(p4.p4)");
      tie(edis,inserted) =   add_edge(1, 3, t1); 
      t1[edis].gmomentum = get_symbol("p5");
      t1[edis].gexpr = 1/get_symbol("(p5.p5)");
      out_vec.push_back(std::make_pair(t1,Tt1));

  
//*************** Topology T2****************
//                   [P3] 
//                .---<-----.
//             3 /           \ 2
//              /\           /\
//          [P1]| *----<----*  | 
//        [Q]___|     [P4]     |___
//              |0           1 |
//              \             /  
//               \_____/____ /
//                     \ 
//                   [P2]
//
//******************************************
      Graph t2;
      tie(edis,inserted) =   add_edge(3, 0, t2); 
      t2[edis].gmomentum = get_symbol("p1");
      t1[edis].gexpr = 1/get_symbol("(p1.p1)");
      tie(edis,inserted) =   add_edge(1, 0, t2);
      t2[edis].gmomentum = get_symbol("p2"); 
      t1[edis].gexpr = 1/get_symbol("(p2.p2)");
      tie(edis,inserted) =   add_edge(2, 3, t2);
      t2[edis].gmomentum = get_symbol("p3"); 
      t1[edis].gexpr = 1/get_symbol("(p3.p3)");
      tie(edis,inserted) =   add_edge(2, 3, t2);
      t2[edis].gmomentum = get_symbol("p4"); 
      t1[edis].gexpr = 1/get_symbol("(p4.p4)");
      tie(edis,inserted) =   add_edge(1, 2, t2);
      t2[edis].gmomentum = get_symbol("p1"); 
      t1[edis].gexpr = 1/get_symbol("(p1.p1)");
      out_vec.push_back(std::make_pair(t2,Tt2));
      break;
      }
    case 3:
      {
        //topology LA
        Graph la;
        tie(edis,inserted) = add_edge(2, 3, la); 
        la[edis].gmomentum = get_symbol("p1");
        la[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(1, 2, la); 
        la[edis].gmomentum = get_symbol("p2");
        la[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(0, 1, la); 
        la[edis].gmomentum = get_symbol("p3");
        la[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(5, 0, la); 
        la[edis].gmomentum = get_symbol("p4");
        la[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(4, 5, la); 
        la[edis].gmomentum = get_symbol("p5");
        la[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(3, 4, la); 
        la[edis].gmomentum = get_symbol("p6");
        la[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(4, 2, la); 
        la[edis].gmomentum = get_symbol("p7");
        la[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(1, 5, la); 
        la[edis].gmomentum = get_symbol("p8");
        la[edis].gexpr = 1/get_symbol("(p8.p8)");

        out_vec.push_back(std::make_pair(la,Tla));

        //topology BE
        Graph be;
        tie(edis,inserted) = add_edge(2, 3, be); 
        be[edis].gmomentum = get_symbol("p1");
        be[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(1, 2, be); 
        be[edis].gmomentum = get_symbol("p2");
        be[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(0, 1, be); 
        be[edis].gmomentum = get_symbol("p3");
        be[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(4, 0, be); 
        be[edis].gmomentum = get_symbol("p4");
        be[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(3, 4, be); 
        be[edis].gmomentum = get_symbol("p5");
        be[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(2, 5, be); 
        be[edis].gmomentum = get_symbol("p6");
        be[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(5, 1, be); 
        be[edis].gmomentum = get_symbol("p7");
        be[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(5, 4, be); 
        be[edis].gmomentum = get_symbol("p8");
        be[edis].gexpr = 1/get_symbol("(p8.p8)");

        out_vec.push_back(std::make_pair(be,Tbe));

          //topology NO
        Graph no;
        tie(edis,inserted) = add_edge(2, 3, no); 
        no[edis].gmomentum = get_symbol("p1");
        no[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(1, 2, no); 
        no[edis].gmomentum = get_symbol("p2");
        no[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(0, 1, no); 
        no[edis].gmomentum = get_symbol("p3");
        no[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(5, 0, no); 
        no[edis].gmomentum = get_symbol("p4");
        no[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(4, 5, no); 
        no[edis].gmomentum = get_symbol("p5");
        no[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(3, 4, no); 
        no[edis].gmomentum = get_symbol("p6");
        no[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(2, 5, no); 
        no[edis].gmomentum = get_symbol("p7");
        no[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(4, 1, no); 
        no[edis].gmomentum = get_symbol("p8");
        no[edis].gexpr = 1/get_symbol("(p8.p8)");

        out_vec.push_back(std::make_pair(no,Tno));

   //topology O1
        Graph o1;
        tie(edis,inserted) = add_edge(1, 2, o1); 
        o1[edis].gmomentum = get_symbol("p1");
        o1[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(0, 1, o1); 
        o1[edis].gmomentum = get_symbol("p2");
        o1[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(3, 0, o1); 
        o1[edis].gmomentum = get_symbol("p3");
        o1[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(2, 3, o1); 
        o1[edis].gmomentum = get_symbol("p4");
        o1[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(3, 4, o1); 
        o1[edis].gmomentum = get_symbol("p5");
        o1[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(4, 5, o1); 
        o1[edis].gmomentum = get_symbol("p6");
        o1[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(4, 5, o1); 
        o1[edis].gmomentum = get_symbol("p7");
        o1[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(5, 1, o1);    // same
        o1[edis].gmomentum = get_symbol("p5"); // momentum
        o1[edis].gexpr = 1/get_symbol("(p5.p5)");  // as in(3,4)

        out_vec.push_back(std::make_pair(o1,To1));

     //topology O2
        Graph o2;
        tie(edis,inserted) = add_edge(1, 2, o2); 
        o2[edis].gmomentum = get_symbol("p1");
        o2[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(0, 1, o2); 
        o2[edis].gmomentum = get_symbol("p2");
        o2[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(5, 0, o2); 
        o2[edis].gmomentum = get_symbol("p3");
        o2[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(4, 5, o2); 
        o2[edis].gmomentum = get_symbol("p4");
        o2[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(5, 1, o2); 
        o2[edis].gmomentum = get_symbol("p5");
        o2[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(3, 4, o2); 
        o2[edis].gmomentum = get_symbol("p6");
        o2[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(3, 4, o2); 
        o2[edis].gmomentum = get_symbol("p7");
        o2[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(2, 3, o2);    // same
        o2[edis].gmomentum = get_symbol("p4"); // momentum
        o2[edis].gexpr = 1/get_symbol("(p4.p4)");  // as in(4,5)

        out_vec.push_back(std::make_pair(o2,To2));



     //topology O4
        Graph o4;
        tie(edis,inserted) = add_edge(2, 3, o4); 
        o4[edis].gmomentum = get_symbol("p1");
        o4[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(1, 2, o4); 
        o4[edis].gmomentum = get_symbol("p2");
        o4[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(5, 1, o4); 
        o4[edis].gmomentum = get_symbol("p3");
        o4[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(3, 5, o4); 
        o4[edis].gmomentum = get_symbol("p4");
        o4[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(5, 2, o4); 
        o4[edis].gmomentum = get_symbol("p5");
        o4[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(3, 4, o4); 
        o4[edis].gmomentum = get_symbol("p6");
        o4[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(0, 4, o4); 
        o4[edis].gmomentum = get_symbol("p7");
        o4[edis].gexpr = 1/get_symbol("(p7.p7)");
        tie(edis,inserted) =   add_edge(0, 1, o4);   // same momentum as (3,4)
        o4[edis].gmomentum = get_symbol("p6"); //
        o4[edis].gexpr = 1/get_symbol("(p6.p6)");  //

        out_vec.push_back(std::make_pair(o4,To4));

     //topology Y1
        Graph y1;
        tie(edis,inserted) = add_edge(2, 3, y1); 
        y1[edis].gmomentum = get_symbol("p1");
        y1[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(4, 3, y1); 
        y1[edis].gmomentum = get_symbol("p2");
        y1[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(1, 2, y1); 
        y1[edis].gmomentum = get_symbol("p3");
        y1[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(1, 2, y1); 
        y1[edis].gmomentum = get_symbol("p4");
        y1[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(5, 4, y1); 
        y1[edis].gmomentum = get_symbol("p5");
        y1[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(5, 4, y1); 
        y1[edis].gmomentum = get_symbol("p6");
        y1[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(0, 1, y1); 
        y1[edis].gmomentum = get_symbol("p1"); // same as (2,3)
        y1[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(0, 5, y1); 
        y1[edis].gmomentum = get_symbol("p2");//same as (4,3)
        y1[edis].gexpr = 1/get_symbol("(p2.p2)");

        out_vec.push_back(std::make_pair(y1,Ty1));

     //topology Y2
        Graph y2;
        tie(edis,inserted) = add_edge(2, 3, y2); 
        y2[edis].gmomentum = get_symbol("p1");
        y2[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(0, 5, y2); 
        y2[edis].gmomentum = get_symbol("p2");
        y2[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(3, 4, y2); 
        y2[edis].gmomentum = get_symbol("p3");
        y2[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(3, 4, y2); 
        y2[edis].gmomentum = get_symbol("p4");
        y2[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(1, 2, y2); 
        y2[edis].gmomentum = get_symbol("p5");
        y2[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(1, 2, y2); 
        y2[edis].gmomentum = get_symbol("p6");
        y2[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(0, 1, y2); 
        y2[edis].gmomentum = get_symbol("p1"); //same as (2,3)
        y2[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(4, 5, y2); 
        y2[edis].gmomentum = get_symbol("p1"); // same
        y2[edis].gexpr = 1/get_symbol("(p1.p1)");

        out_vec.push_back(std::make_pair(y2,Ty2));

     //topology Y3
        Graph y3;
        tie(edis,inserted) = add_edge(4, 5, y3); 
        y3[edis].gmomentum = get_symbol("p1");
        y3[edis].gexpr = 1/get_symbol("(p1.p1)");
        tie(edis,inserted) =   add_edge(0, 5, y3); 
        y3[edis].gmomentum = get_symbol("p2");
        y3[edis].gexpr = 1/get_symbol("(p2.p2)");
        tie(edis,inserted) =   add_edge(3, 4, y3); 
        y3[edis].gmomentum = get_symbol("p3");
        y3[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(1, 4, y3); 
        y3[edis].gmomentum = get_symbol("p4");
        y3[edis].gexpr = 1/get_symbol("(p4.p4)");
        tie(edis,inserted) =   add_edge(2, 3, y3); 
        y3[edis].gmomentum = get_symbol("p5");
        y3[edis].gexpr = 1/get_symbol("(p5.p5)");
        tie(edis,inserted) =   add_edge(2, 3, y3); 
        y3[edis].gmomentum = get_symbol("p6");
        y3[edis].gexpr = 1/get_symbol("(p6.p6)");
        tie(edis,inserted) =   add_edge(1, 2, y3); 
        y3[edis].gmomentum = get_symbol("p3"); // same as (3,4)
        y3[edis].gexpr = 1/get_symbol("(p3.p3)");
        tie(edis,inserted) =   add_edge(0, 1, y3); 
        y3[edis].gmomentum = get_symbol("p1");//same as (4,5)
        y3[edis].gexpr = 1/get_symbol("(p1.p1)");

        out_vec.push_back(std::make_pair(y3,Ty3));






      }
      break;
    }
      return out_vec;
    }


/*
std::vector<Graph> t2_full()
{
  Graph t1;
  add_edge(0, 1, t1); 
  add_edge(1, 3, t1); 
  //  add_edge(1, 3, t1); 
  add_edge(3, 2, t1); 
  add_edge(2, 0, t1); 
  add_edge(2, 1, t1); 
  return t1;
}
*/
 // shrink 1 edge
Graph shrink_edge(Graph::edge_descriptor edis_in,Graph G)
{
  Graph new_g;
  Graph::edge_iterator ei,ei_end;
  Graph::vertex_descriptor u,v,src,tar;
  Graph::edge_descriptor edis;
  // shrink to source vertex
  u = source(edis_in,G);
  v = target(edis_in,G);
  bool inserted;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    {
      if(v==source(*ei,G))
        {
          src = u;
          tar = target(*ei,G);
        }
      else if(v==target(*ei,G))
        {
          tar = u;
          src = source(*ei,G);
        }
      else
        {
          src = source(*ei,G);
          tar = target(*ei,G);
        }
      
      tie(edis, inserted) = add_edge(src,tar,new_g);
      new_g[edis].gexpr = G[*ei].gexpr;
      new_g[edis].gmomentum = G[*ei].gmomentum;
    }
  return new_g;
}
// shrink list of edges
Graph shrink_edge_list(std::vector<Graph::edge_descriptor> edis_vec,Graph& G)
{
  Graph new_g;
  cout << " Edges to remove "<<edis_vec.size()<<"  nu ed start:  "<<num_edges(G)<<endl;
  std::vector<std::pair<Graph::vertex_descriptor,Graph::vertex_descriptor> > shrink_pair_vec;
  BOOST_FOREACH(Graph::edge_descriptor ed,edis_vec)
    {
      shrink_pair_vec.push_back(make_pair(source(ed,G),target(ed,G))); //important step
      if(find(edis_vec.begin(),edis_vec.end(),ed)!=edis_vec.end())
        cout<<ed<<" found"<<endl;
      else 
        cout<<ed<<" NOT found"<<endl;
    }



  print_graph(G);
  typedef  std::map<Graph::edge_descriptor,std::pair<Graph::vertex_descriptor,Graph::vertex_descriptor> > edge_map;
  edge_map emap;
  Graph::edge_iterator ei,ei_end;
  Graph::vertex_descriptor u,v,src,tar;
  Graph::edge_descriptor edis;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    {
      Graph::edge_descriptor ert = *ei;
      if(find(edis_vec.begin(),edis_vec.end(),ert)!=edis_vec.end())
        cout<<ert<<" found"<<endl;
      else 
        {
          cout<<ert<<" NOT found"<<endl;
          emap[*ei] = std::make_pair(source(*ei,G),target(*ei,G));
        }
    }


  // shrinking second to first
  for(std::vector<std::pair<Graph::vertex_descriptor,Graph::vertex_descriptor> >::iterator it = shrink_pair_vec.begin();it!=shrink_pair_vec.end();++it)
    {
      for(std::vector<std::pair<Graph::vertex_descriptor,Graph::vertex_descriptor> >::iterator it_v = boost::next(it);it_v!=shrink_pair_vec.end();++it_v) // modyfying shrink list
          {
            if(it_v->second == it->second)it_v->second = it->first;
            if(it_v->first == it->second)it_v->first = it->first;
          }
      for(edge_map::iterator mit = emap.begin();mit!=emap.end();++mit)
        {
          if(mit->second.second == it->second)mit->second.second = it->first;
          if(mit->second.first == it->second)mit->second.first = it->first;
        }
    }
  bool inserted;
  for(edge_map::iterator mit = emap.begin();mit!=emap.end();++mit)
    {
      cout<<">> to ins: "<<mit->second.first<<"  ,  "<<mit->second.second<<endl;
      tie(edis, inserted) = add_edge(mit->second.first,mit->second.second,new_g);
      new_g[edis].gexpr = G[mit->first].gexpr;
      new_g[edis].gmomentum = G[mit->first].gmomentum;
      
    }

  print_graph(new_g);

  /*  
// shrink to source vertex
  u = source(edis_in,G);
  v = target(edis_in,G);
  bool inserted;
  for(tie(ei,ei_end)=edges(G);ei!=ei_end;++ei)
    {
      if(v==source(*ei,G))
        {
          src = u;
          tar = target(*ei,G);
        }
      else if(v==target(*ei,G))
        {
          tar = u;
          src = source(*ei,G);
        }
      else
        {
          src = source(*ei,G);
          tar = target(*ei,G);
        }
      
      tie(edis, inserted) = add_edge(src,tar,new_g);
      new_g[edis].gexpr = G[*ei].gexpr;
    }

  */
  return new_g;
}



// get vector of graph with n_v vertices and n_e edges
Graph_topo get_topo_v_e(int n_v,int n_e)
{
  try
    {
      Graph_topo out_vec;
      Graph_topo in_vec(get_root_topos(n_e-n_v+1)); // loops in 1PI diagram
      cout<<"Get topo started"<<endl;
      BOOST_FOREACH(Graph_topo_pair graph_topo_cntr,in_vec)
        {
          Graph root_g(graph_topo_cntr.first);
          int n = num_edges(root_g);
          int r = num_edges(root_g) - n_e;
          if((n_v == num_vertices(root_g) ) && (n_e == num_edges(root_g) ))
            out_vec.push_back(graph_topo_cntr);
          else
            {
              Graph::edge_iterator ei,ei_end;
              std::vector<int> v_int(n);
              for (int i = 0; i < n; ++i) { v_int[i] = i; }
              std::map<int,Graph::edge_descriptor> edis_map;
              int cntr = 0;
              for(tie(ei,ei_end)=edges(root_g);ei!=ei_end;++ei)  
                {
                  edis_map[cntr] = *ei;
                  cntr++;
                }
  
              do{
                std::vector<Graph::edge_descriptor> ed_vec;
                for (int j = 0; j < r; ++j) 
                  { 
                    std::cout << " >>edge  " <<edis_map[ v_int[j]]; 
                    ed_vec.push_back(edis_map[ v_int[j]]);
                  }
                cout<<endl<<endl;
                Graph tm_gr = shrink_edge_list(ed_vec,root_g);
                if((num_vertices(tm_gr) == n_v) && (num_edges(tm_gr) == n_e)) // only graphs with n_e and n_v are accepted
                  out_vec.push_back(std::make_pair(tm_gr,graph_topo_cntr.second));
    
              }while (next_combination(v_int.begin(), v_int.begin() + r, v_int.end()));
            }
        }
  return out_vec;
}
  catch(std::exception& p)
    {
      //      cout<<"*********************************************"<<endl;
      // cout<<">>>ERROR:  "<<setw(40)<<std::internal<<what<<endl;
      //cout<<"*********************************************"<<endl;
      throw std::logic_error(std::string("In function \"get_topo_v_e\":\n |___> ")+p.what());
    }
}



